function newFunctionValues = linearParabolicPDE_timeStep(spaceGrid,functionValues,dt,...
  diffusion,advection,discounting,inhomogeneity,boundaryLeft,boundaryRight)
%SOLVEKFE_OU Performs a single forward time step of the following PDE
%   v_t = diffusion*v_xx + advection*v_x + discounting*v + inhomogeneity
% solver uses an implicit Euler scheme and upwind differencing for first
% derivatives

  nNodesSpace = numel(spaceGrid);
  method = 1; % 1 is implicit, code below can also handle explicit and mixed
  
  % handling of differentiation formulas in code:
  % variable naming: cD1l, cD1m, cD1r are coefficients in differentiation
  % formula
  %   f'(x_i) ~ cD1l*f(x_{i-1}) + cD1m*f(x_i) + cD1r*f(x_{i+1})
  % --------------------------------
  % as these coefficients depend on the position in the grid, cF1l, cF1m,
  % cD1r are vectors indexed by the grid point index i
  % --------------------------------
  
  cD1ll = [NaN; - 1./(spaceGrid(2:end-1)-spaceGrid(1:end-2)); NaN];
  cD1ml = [NaN; 1./(spaceGrid(2:end-1)-spaceGrid(1:end-2)); NaN];
  cD1rl = [NaN; zeros(nNodesSpace-2,1); NaN];
  cD1lr = [NaN; zeros(nNodesSpace-2,1); NaN];
  cD1mr = [NaN; - 1./(spaceGrid(3:end)-spaceGrid(2:end-1)); NaN];
  cD1rr = [NaN; 1./(spaceGrid(3:end)-spaceGrid(2:end-1)); NaN];
  posCoeff = advection>=0;
  negCoeff = advection<0;
  cD1l = posCoeff.*cD1lr + negCoeff.*cD1ll;
  cD1m = posCoeff.*cD1mr + negCoeff.*cD1ml;
  cD1r = posCoeff.*cD1rr + negCoeff.*cD1rl;
  
  % second order differentiation formula
  cD2l = [NaN; 2./(spaceGrid(3:end)-spaceGrid(1:end-2))./(spaceGrid(2:end-1)-spaceGrid(1:end-2)); NaN];
  cD2m = [NaN; -2./(spaceGrid(3:end)-spaceGrid(1:end-2)).*(1./(spaceGrid(3:end)-spaceGrid(2:end-1)) + 1./(spaceGrid(2:end-1)-spaceGrid(1:end-2))); NaN];
  cD2r = [NaN; 2./(spaceGrid(3:end)-spaceGrid(1:end-2))./(spaceGrid(3:end)-spaceGrid(2:end-1)); NaN];

  newFunctionValues = NaN(nNodesSpace,1);
  newFunctionValues(1) = boundaryLeft;
  newFunctionValues(end) = boundaryRight;

  % prepare equation matrix
  % equation system for implicit method is
  %   (I-dt*equMatrix)*v = rhs
  % where dt is the time increment and rhs is a vector defined in the 
  % iteration
  % ------------------------------
  % the more general equation for the implicit/explicit/mixed method is
  %   (I-alpha*dt*equMatrix)*p_{n+1} = (I+(1-alpha)*dt*equMatrix)*p_{n} + dt*inhomogeneity
  % where alpha is the method parameter in [0,1] and p_{n}, p_{n+1} are
  % vectors of the solution function at the last and new time.
  % (at the boundary rows, the equation must be slightly corrected for the
  % presence of boundary conditions)
  equMatrix = spalloc(nNodesSpace-2,nNodesSpace-2,3*(nNodesSpace-2)-2);
  equMatrixD = discounting(2:end-1) + advection(2:end-1).*cD1m(2:end-1) + diffusion(2:end-1).*cD2m(2:end-1); %diagonal part of matrix
  equMatrixUD = advection(2:end-1).*cD1r(2:end-1) + diffusion(2:end-1).*cD2r(2:end-1); %upper diagonal
  equMatrixLD = advection(2:end-1).*cD1l(2:end-1) + diffusion(2:end-1).*cD2l(2:end-1); %lower diagonal
  equMatrix = spdiags([equMatrixD,[equMatrixLD(2:end);NaN],[NaN;equMatrixUD(1:end-1)]],[0,-1,1],equMatrix);
  % do iteration (implicit only for now)
  equRHS = (speye(nNodesSpace-2)+(1-method)*dt*equMatrix)*functionValues(2:end-1) + dt*inhomogeneity(2:end-1);
  equRHS(1) = equRHS(1) + method*dt*equMatrixLD(1)*newFunctionValues(1); %correct leftmost term by knowledge about left boundary
  equRHS(end) = equRHS(end) + method*dt*equMatrixUD(end)*newFunctionValues(end); %correct rightmost term by knowledge about right boundary
  % solve equation system for implicit portion
  newFunctionValues(2:end-1) = (speye(nNodesSpace-2)-method*dt*equMatrix)\equRHS;
end

